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Abstract 

Performance  equations  that  describe  the  dependence  of  cathode  potential  on  current  density  for  polymer  electrolyte  fuel  cells  (PEFCs) 
are  developed  based  on  a  mechanistic  approach.  The  equations  take  into  account,  in  detail,  potential  losses  caused  by:  (i)  electric  resistance 
and  gas  transport  limitations  in  the  gas  diffuser;  (ii)  limitations  of  oxygen  diffusion,  proton  migration  and  electron  conduction  in  the  catalyst 
layer;  (iii)  oxygen  reduction  within  the  catalyst  layer.  Derivation  of  the  equations  is  initiated  with  the  formulation  of  a  one-dimensional 
model  and  followed  by  the  incorporation  of  appropriate  profiles  for  oxygen  concentration,  ionomer  potential  and  catalyst  potential.  The 
final  forms  of  the  equations  are  obtained  by  lumping  the  oxygen  reduction  reaction  at  a  reaction  center  of  the  catalyst  layer.  Since  the 
equations  are  derived  from  a  mechanistic  model,  all  parameters  appearing  in  the  equations  are  endowed  with  a  precise  physical  meaning. 
In  addition,  potential  losses  caused  by  various  sources  can  be  clearly  quantified.  Excellent  agreement  is  found  between  the  results  obtained 
from  the  equations  and  from  the  one-dimensional  model  over  an  extensive  range  of  the  values  of  model  parameters.  This  indicates  that 
the  present  equations  can  be  employed  to  replace  the  one-dimensional  model  as  an  efficient  and  convenient  predictive  tool  for  cathode 
performance  with  greatly  reduced  computational  efforts  while  keeping  the  same  level  of  accuracy. 
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1.  Introduction 

Due  to  their  promising  application  as  a  power  genera¬ 
tion  device  in  electric  vehicles,  polymer  electrolyte  fuel  cells 
(PEFCs)  have  received  increasing  attention  during  the  past 
decade.  Parallel  to  the  enhancements  in  the  efficiency  and 
stability  of  PEFCs  by  the  development  of  improved  compo¬ 
nents  and  fabrication  techniques,  progress  in  mathematical 
modelling  has  provided  a  deeper  insight  into  the  physical 
and  chemical  processes  as  an  aid  to  better  design  and  ap¬ 
propriate  operation. 

Mathematical  models  with  varying  degrees  of  sophisti¬ 
cation  in  transport  phenomena  and  spatial  dimensionality 
have  been  developed.  Rigorous  formulation  of  a  complete 
cell  with  the  membrane-electrode-assembly  and  the  gas  flow 
channels  leads  to  a  three-dimensional  model  [1,2].  In  ad¬ 
dition  to  oxygen  transport  normal  to  the  gas-diffuser  face, 
different  considerations  arise  with  a  two-dimensional  de¬ 
sign,  such  as  variation  in  oxygen  concentration  along  the 
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gas  channels  [3-6],  inaccessibility  on  the  face  of  the  gas 
diffuser  to  oxygen  gases  due  to  the  rib  structure  of  carbon 
plates  [7-9],  or  interdigitated  flow  fields  [10].  With  differ¬ 
ent  assumptions,  wide  variations  can  also  be  found  in  the 
one-dimensional  PEFC  models  [1 1-18].  Without  taking  spa¬ 
tial  gradients  of  concentration  and  potential  into  account, 
PEFC  models  in  the  form  of  algebraic  equations  have  also 
been  proposed  [19-24]. 

Based  on  the  methodologies  in  formulation,  the  PEFC 
models  can  be  considered  as  mechanistic,  empirical,  or 
semi-empirical  in  nature.  A  mechanistic  model  is  derived 
from  phenomenological  mass  transport  and  conservation 
equations.  Thus,  the  complexity  of  the  models  varies  with 
the  mechanisms  taken  into  account.  On  the  other  hand,  an 
empirical  model  is  typically  one  equation  composed  of  a 
few  terms  that  account  for  the  observed  characteristics  of 
experimental  data.  A  semi-empirical  model  is  featured  by 
the  characteristics  of  both  approaches.  A  complex  mecha¬ 
nistic  model  is  formulated  in  a  more  involved  way  and  thus 
it  is  able  not  only  to  exhibit  a  great  predictive  power  but  to 
provide  a  better  interpretation  of  experimental  data  as  well. 
In  addition  to  the  considerable  modelling  and  computational 
efforts  that  are  required,  the  model  parameters  governing 
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Nomenclature 

a  parameter  in  Eq.  (24) 

A  effective  platinum  surface  area  per  unit 
volume  (cm2  cm-3) 

co2  oxygen  concentration  (mol  cm-3) 

co2,ref  reference  oxygen  concentration 

(mol  cm-3),  defined  as  P/Ho2 
d  gas-diffuser  thickness  (mm) 

D1-;  effective  binary  diffusion  coefficient 

for  i  and  j  species  (cm2  s-1) 
effective  diffusivity  of  dissolved  oxygen  in 
the  catalyst  layer  (cmz  s  ) 

F  Faraday’s  constant  (96,487  C  per 

equivalent) 

Ho2  Henry’s  constant  for  oxygen 

(atm  cm3  mol-1) 

i  protonic  current  density  (A  cm-2) 

/e  electronic  current  density  (A  cm-2) 

4, ref  exchange  current  density  at  reference 
condition  (A  cm-2) 

I  cathode  current  density  (A  cm-2) 

IQ  reference  current  density  (A  cm-2), 

defined  by  Eq.  (Ilf) 
effective  electronic  conductivity  for 
gas  diffuser  (mho  cm-1) 
k effective  protonic  conductivity  for 
the  ionomer  phase  (mho  cm-1) 
kf^  effective  electronic  conductivity  for  solid 

catalysts  (mho  cm-1) 

n  number  of  electrons  transferred  in  reaction 

Ni  molar  flux  of  species  i  (mol  cm-2) 

P  total  pressure  (atm) 

R  the  universal  gas  constant 

(8.314  Jmol-1K-1) 

T  cathode  temperature  (K) 

V  potential  of  catalysts  (V) 

V*  catalyst  potential  at  reaction  center  (V) 

yc  cathode  potential  (V) 

Vq  equilibrium  potential  (V) 

w  dimensionless  position  where  oxygen 

depletion  occurs 

Xi  molar  fraction  of  species  i 

x02  average  value  of  xq2  over  catalyst  layer 
Vq2  xo2  at  catalyst-layerlgas-diffuser  interface 

z  co-ordinate  perpendicular  to  face  of 

the  gas  diffuser  (pan) 

Greek  letters 

oti  electrode  transfer  coefficient 

Pi  parameter,  defined  by  Eq.  (lib)  (mho cm-2) 

P2  parameter,  defined  by  Eq.  (lie)  (mho cm-2) 

Pi  parameter,  defined  by  Eq.  (28)  (mho  cm-2) 

P\  parameter,  defined  by  Eq.  (33)  (A-1  cm2) 

y  kinetic  coefficient  taken  as  unity 


8  the  catalyst-layer  thickness  (pan) 

f  dimensionless  co-ordinate,  defined 

by  Eq.  (lid) 

0  ionomer  potential  (V) 

0*  ionomer  potential  at  the  reaction  center  (V) 
cp  parameter,  defined  by  Eq.  (11a) 


the  various  sources  of  polarization  within  the  PEFCs  are 
related  in  a  complex  nature  that  cannot  be  easily  captured. 
In  other  respects,  empirical  models  are  characterized  by  a 
simple  algebraic  expression  that  relates  cell  potential  and 
current  density  [19-21].  Due  to  the  lack  of  physical  origin 
and  significance  in  the  model  parameters,  however,  they  are 
abundant  with  fitting  coefficient  which  thus  eliminates  the 
possibility  of  using  them  as  predictive  tools.  In  order  to  cir¬ 
cumvent  the  drawbacks  of  empirical  models,  semi-empirical 
formulation  of  a  performance  equation  has  been  proposed 
to  reduce  the  number  of  fitting  coefficients,  in  which  both 
mechanistic  and  empirical  features  are  present  in  the  model 
parameters  [22-25]. 

In  the  present  work,  an  attempt  is  made  to  derive  perfor¬ 
mance  equations  with  all  mechanistic  derived  coefficients. 
In  order  to  do  this,  attention  is  first  paid  to  the  chemical 
and  transport  phenomena  that  occur  in  the  cathode  of  a 
PEFC  and  formulate  a  one-dimensional  mechanistic  model 
that  takes  into  account  spatial  variations  of  concentration 
and  potential  in  the  catalyst  layer  and  the  gas  diffuser.  By 
approximating  the  concentration  and  potential  profiles  with 
appropriate  shape  functions  and  lumping  the  oxygen  reduc¬ 
tion  reaction  at  the  reaction  center,  it  is  possible  to  reduce 
the  one-dimensional  model  to  a  single  performance  equa¬ 
tion  that  presents  the  relationship  between  cathode  potential 
and  current  density.  Because  the  performance  equations  are 
derived  from  a  mechanistic  model,  they  can  be  employed 
as  a  predictive  tool  while  characterized  by  simple  algebraic 
expression  just  like  empirical  model  equations.  In  addition, 
individual  contributive  terms  to  the  overall  cathode  polariza¬ 
tion  can  be  clearly  identified.  Computational  results  also  re¬ 
veal  that  excellent  agreement  exists  between  the  polarization 
curves  (calculated  based  on  the  one-dimensional  model)  and 
the  developed  performance  equations.  This  indicates  that  the 
derived  performance  equations  can  be  used  as  a  replacement 
for  the  one-dimensional  model  without  introducing  any  ap¬ 
preciable  errors. 

2.  One-dimensional  model 

The  mathematical  model  for  the  PEFC  cathodes  consid¬ 
ered  in  the  present  work  is  a  one-dimensional,  steady-state 
and  isothermal  model.  The  cathodes  are  composed  of  two 
major  components,  namely,  a  catalyst  layer  and  a  gas  dif¬ 
fuser.  In  the  former,  it  is  assumed  that  catalyst  particles,  usu¬ 
ally  in  the  form  of  platinum  supported  on  carbon  black  or  of 
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unsupported  platinum  black,  are  homogeneously  mixed  with 
a  proton-conductive  ionomer.  The  gas  diffuser,  which  is  in 
intimate  contact  with  the  catalyst  layer,  is  made  of  carbon  pa¬ 
per  or  carbon  cloth  with  hydrophobic  treatment.  During  op¬ 
eration,  oxygen  gases,  mostly  coming  from  a  humidified  air 
stream,  penetrate  the  gas  diffuser  to  reach  its  interface  with 
the  catalyst  layer,  and  then  dissolve  into  the  ionomer  phase. 
Since  the  ionomer  consists  of  both  hydrophobic  branches 
and  hydrophilic  groups,  there  might  also  exist  some  void 
spaces  in  the  catalyst  layer.  Thus,  some  oxygen  gases  will 
diffuse  further  into  the  inner  region  of  the  catalyst  layer  via 
gas  channels,  and  dissolve  into  the  ionomer  phase  thereby. 
In  the  ionomer  phase,  the  dissolved  oxygen  will  continue  to 
travel  to  the  surface  of  the  platinum  particles,  where  an  elec¬ 
trochemical  reaction  with  the  protons  from  the  anode  occurs 
and  can  be  described  by 

4H+  +  4e“  +  02  2H20 


The  Bulter-Volmer  equation  for  the  rate  of  the  above  elec¬ 
trochemical  reaction  gives 


V-i  —  A  4,  ref 


co2 

c02,ref 


—  exp 


where  oti  is  the  electrode  transfer  coefficient;  i  the  protonic 
current  density;  A  the  effective  platinum  surface  area  per  unit 
volume;  4, ref  is  the  exchange  current  density  at  the  reference 
condition;  y  a  kinetic  coefficient  and  taken  as  unity;  T  the 
temperature  of  the  catalyst  layer;  R  the  universal  gas  con¬ 
stant;  co2  the  dissolved  oxygen  concentration;  n  the  number 
of  electrons  transferred  in  the  reaction;  V0  the  equilibrium 
potential;  V  the  catalyst  potential;  and  0  the  ionomer  poten¬ 
tial. 

The  mass-balance  equation  of  oxygen  in  the  catalyst  layer 
can  be  written  as 

<y2c°2  +  ^  =  0  (2) 

in  which  is  the  effective  diffusivity  of  dissolved  oxy¬ 
gen  in  the  catalyst  layer.  As  stated  previously,  oxygen  may 
diffuse  in  the  catalyst  layer  via  different  paths,  i.e.  gas  pores 
or  hydrated  ionomer.  In  fact,  the  value  of  may  vary 
significantly  due  to  the  various  degrees  of  involvement  of 
these  two  different  diffusion  routes.  In  the  present  study,  the 
details  of  the  oxygen  diffusion  mechanism  are  not  consid¬ 
ered.  Instead,  is  treated  as  an  adjustable  parameter  so 
as  to  investigate  its  influence  on  cathode  performance  and 
to  explore  the  discrepancies  between  the  performance  equa¬ 
tions  and  the  one-dimensional  model  under  different  rates 
of  oxygen  transport.  The  ohmic  potential  losses  in  the  cat¬ 
alyst  layer  are  attributed  to  both  proton  migration  through 
the  ionomer  phase  and  electron  conduction  in  the  solid  cata¬ 
lysts.  Applying  Ohm’s  law  to  the  motion  of  both  the  proton 


and  the  electrons  gives 
dV  4 

—  =  -4?  (3) 

dz  kf1 

and 

d  0  i 

—  = - -  (4) 

dz  k^f 

where  kff  and  k^  denote  the  effective  electric  conductivity 
for  the  solid  catalysts  and  the  effective  protonic  conductiv¬ 
ity  for  the  ionomer  phase,  respectively,  and  4  represents  the 
electronic  current  density.  From  electroneutrality,  one  ob¬ 
tains 

V  •  i  +  V  •  /e  =  0  (5) 


As  a  humidified  air  stream  is  used  as  the  cathode  feed,  three 
gaseous  species,  viz.  oxygen,  nitrogen  and  water  vapor,  exist 
in  the  diffusion  layer.  The  Stefan-Maxwell  equations  for 
multi-component  diffusion  are  usually  a  starting  point  for 
the  description  of  gas  transport  in  the  gas  diffuser,  and  take 
the  form 


Vxi  = 


RT 


n 

one ff 

7=1  i-J 


( XiNj  —  XjNi),  i  =  N2,  w,  02 


where  is  an  effective  binary  diffusion  coefficient  in 
the  porous  medium  for  i  and  j  species  and  Nj  the  molar 
flux  of  species  i  with  jc*  its  molar  fraction.  It  is  assumed 
that  a  phase  equilibrium  between  the  water  vapor  and  its 
condensed  phase  is  achieved  at  each  point  of  the  gas  diffuser. 
Also,  the  molar  flux  of  nitrogen  is  taken  to  be  zero  due  to 
its  inertness.  Thus,  the  Stefan-Maxwell  equations  can  be 
condensed  into  one  single  equation  [13],  i.e. 


xo2)no2 


x 


1 


Deff 

L  ^n2-o2 


+ 


w 


*02£>N2-w  +  <7 


*o2><-wJ 

(7) 


In  order  to  facilitate  computations,  the  above  equations  were 
rearranged  and  written  in  partially  dimensionless  form.  In 
the  catalyst  layer  (0  <  f  <  1) 


d2xp2 

d£2 


V 


exp 


exp 


acnF(Vo  -  V  +  0) 

RT 

—aanF(Vo  —  V  +  <j>) 
RT 


x0,  =  0 


(i\  d2V  d2xp2  _  o 


4>  dc 


df- 


(8) 

(9) 


d2x0z  _  fe  d20 
d£2  I0  d£2 


(10) 
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where 


<P  = 


^o,ref  Ho2$ 
nFPDf 


Pi  = 


Pl  = 


k 


eff 


5 

ceff 

Km 


Z 

^=8 


Xq2 


Hq2co2 


Io  = 


nFPDf 
Ho,  8 


(Ha) 

(lib) 

(11c) 

(lid) 

(lie) 

(Ilf) 


where  8  denotes  the  catalyst-layer  thickness  and  Ho2  the 
Henry’s  constant  for  gaseous  oxygen  and  its  dissolved  form 
in  the  ionomer  phase  at  the  cathode  temperature.  It  should 
be  noted  that  xo2  denotes  the  molar  fraction  of  oxygen  in 
the  gas  phase  of  the  diffuser  but  stands  for  the  dimension¬ 
less  concentration  of  dissolved  oxygen  in  the  ionomer  phase 
within  the  catalyst  layer  as  defined  by  Eq.  (lie).  The  same 
notation  is  used  in  these  two  domains  because  by  such  def¬ 
initions  xo2  profiles  are  continuous  across  their  boundary. 

In  the  gas  diffuser  (1  <  K  <  b),  two  more  equations  are 
required  in  addition  to  Eq.  (7),  viz. 


d/Vo2 

d? 


d2V 

df2 


(12) 

(13) 


where  V  is  the  potential  of  the  carbon  paper  (or  cloth).  This 
is  the  same  symbol  as  that  used  for  the  potential  of  solid 
catalysts  in  the  catalyst  layer  since  these  two  quantities  are 
identical  at  the  boundary  between  these  two  regions.  At  the 
face  of  the  gas  diffuser  (£  =  b ),  two  more  boundary  condi¬ 
tions  are  required  apart  from  Eq.  (12)  i.e. 

xo2  =  xo2  (14) 

V  =  Vc  (15) 


where  indicates  the  oxygen  molar  fraction  in  the  bulk 
flow  and  Vc  is  the  cathode  potential.  At  the  gas-diffuserl 
catalyst-layer  interface,  six  boundary  conditions  are  needed 
since  there  are  three  state  variables  in  each  domain.  First, 
Eq.  (12)  is  also  applied.  In  addition,  continuities  of  V  and 
xo2  require 

xo2  (catalyst  layer)  =  xo2  (gas  diffuser)  (16) 

V  (catalyst  layer)  =  V  (gas  diffuser)  (17) 


Equivalence  of  oxygen  molar  flux  and  electric  current  den¬ 
sity  on  both  sides  of  the  boundary  leads  to 


—  (catalyst  layer)  =  k~f  —  (gas  diffuser) 


(18) 


Io  dvp2 
nF  df 


(catalyst  layer)  =  —  Nq2 


(19) 


Because  protons  are  not  able  to  penetrate  the  catalyst-layer  I 
gas-diffuser  interface,  zero  protonic  current  is  expected  here, 
which  gives 


(20) 


At  the  membrane  I  cataly  st-layer  interface  (£  =  0),  the  bound¬ 
ary  conditions  are  formulated  as 


0  =  0 


(21) 


(22) 


(23) 


The  mathematical  model  described  above  was  formulated 
following  the  lines  of  Bernardi  and  Verbrugge  [13]  as  well 
as  Springer  et  al.  [17].  Compared  with  the  work  of  Bernardi 
and  Verbrugge  [13],  a  membrane  is  not  included  and  flow 
of  liquid  water  in  the  catalyst  layer  is  not  taken  into  account 
in  the  present  model  for  the  sake  of  simplification.  Finite 
values  of  electric  conductivity  for  the  solid  catalysts  and  the 
gas  diffuser  are  considered  to  account  for  ohmic  potential 
loss  due  to  electronic  movement,  which  are  neglected  in  the 
work  of  Springer  et  al.  [17]. 

The  model  equations  derived  above  form  a  three-point 
boundary  value  problem  with  three  state  variables  in  each 
domain.  The  method  of  collocation  on  finite  elements  based 
on  cubic  B- spline  interpolation  was  employed  for  the  nu¬ 
merical  solutions  of  the  model  equations  [26] .  First,  the  do¬ 
mains  of  the  gas  diffuser  and  the  catalyst  layer  were  divided 
into  several  sub-intervals,  respectively.  Because  of  the  mod¬ 
erate  non-linearity  of  the  Stefan-Maxwell  equations  and  the 
mild  steepness  for  the  oxygen  concentration  profiles  devel¬ 
oped  in  the  gas  diffuser,  only  three  sub-intervals  with  equal 
sizes  were  arranged  in  this  domain.  Twelve  sub-intervals  of 
non-equal  sizes  were  employed  in  the  catalyst-layer  domain 
to  account  for  high  non-linearity  of  the  electrochemical 
rate  expression  for  the  Bulter-Volmer  equation  and  soaring 
oxygen  concentration  gradients  developed  at  high  current 
densities.  In  the  catalyst-layer  domain,  the  sizes  of  the 
sub-intervals  were  decreased  along  the  ^-direction  by  a 
common  ratio  of  0.37.  Continuities  of  the  first  derivatives  of 
all  state  variables  were  required  at  each  breakpoint  point  of 
sub-intervals  except  the  gas-diffuserlcatalyst-layer  interface. 
Two  points  were  collocated  in  each  sub-interval,  which 
thus  transforms  the  model  equations  into  a  set  of  algebraic 
equations.  Newton’s  method  was  employed  to  obtain  the 
numerical  solutions  of  the  resulting  algebraic  equations. 
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3.  Approximate  model  I  for  performance  equations  where 


It  is  assumed  that  the  dimensionless  oxygen  concentration 
profile  in  the  catalyst  layer  can  be  described  by  a  parabolic 
function,  which  has  the  form 

xo2  =  at;2  +  xSQ2-a  (24) 


where  Vqo  denotes  the  value  of  xo2  at  the  catalyst-layerlgas- 
diffuser  interface.  Eq.  (24)  naturally  satisfies  the  bound¬ 
ary  conditions  at  f  =  0  and  1.  By  relating  the  current 
density  with  the  flux  of  dissolved  oxygen  at  the  catalyst- 
layerlgas-diffuser  interface,  a  can  be  expressed  as 


a  = 


I 

2/o 


(25) 


RTd 

nFP 


1 


D 


eff 

N2  O2 


w 


Xy,)D 


eff 

N2— w 


(33) 


After  evaluating  V*,  0*  and  xo2in  Eq.  (30)  with  Eq.  (26)  and 
(27)  and  then  taking  logarithms  on  both  sides,  one  arrives  at 


/  (3642/0-ll/)(1242/0-57) 

”  12801  (3/o42-/)2 

I  37  (442/o-/)(3642/0-137) 

ft  "  128/12  (34270  -  7)2 


RT 

acnF 


In 


7 

<P(xS02I0  ~  //3) 


(34) 


where  7  stands  for  the  cathode  current  density.  The  profiles 
of  ionomer  potential  and  catalyst  potential  within  the  cata¬ 
lyst  layer  can  be  easily  derived  from  Eqs.  (9),  (10)  and  (24) 
together  with  their  associated  boundary  conditions,  which 
are 


+  Vc  + 


ft 


(26) 


(27) 


where 


(28) 


k^f  represents  the  effective  electric  conductivity  of  the  gas 
diffuser  and  d  its  thickness.  Furthermore,  the  anodic  reac¬ 
tion  rate  is  neglected  and  the  rate  of  electrochemical  reac¬ 
tion  occurring  in  the  catalyst  layer  is  lumped  at  the  reaction 
center,  which  is  defined  by 


/pfooz  d£  =  3(4xq2/0  —  7) 
fx02  d?  _8(342/o-/) 


The  overall  potential  loss  in  a  PEFC  cathode  can  be 
viewed  as  a  summation  of  diffusion  overpotential  and  ohmic 
overpotential  for  electron  conduction  in  the  gas  diffuser,  dif¬ 
fusion  overpotential  and  ohmic  potential  losses  for  proton 
migration  and  electron  conduction  in  the  catalyst  layer,  and 
activation  overpotential  for  the  electrochemical  reaction.  It 
is  noted  that  each  individual  contribution  in  the  cathode 
potential  loss  can  be  clearly  quantified  in  Eq.  (34),  which 
are 


diffusion  overpotential  in  gas  diffuser 


Ole nF  yx02  J 

(35) 


activation  overpotential  for  electrochemical  reaction 


(36) 


diffusion  overpotential  in  catalyst  layer 

=  RT  in  (  \ 

nF  n  l  Vq270  -  7/3  J 


(37) 


Based  on  the  two  approximations  stated  above,  one  can  have 


7  —  Aiojef8 


exp 


acnF(V0~  E*+0*) 
RT 


xq2 


(30) 


where  E*  and  0* represent  the  catalyst  potential  and  the 
ionomer  potential  at  the  reaction  center,  respectively.  xo2 
stands  for  the  average  value  of  xo2  over  the  catalyst  layer, 
which  is  calculated  by 


1 

xo2  d;  =  xs02 


I 

37o 


(31) 


By  following  the  approach  of  Pisani  et  al.  [25],  that  is 
to  approximate  7)^q2  with  7)^2_w,  Eq.  (7)  can  be  directly 
solved  and  Vq2  can  be  expressed  as 


xo2  =  xo2  ~  xn2[^pW)  -  1] 


(32) 


ohmic  loss  (electronic)  in  catalyst  layer 
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ohmic  loss  (protonic)  in  catalyst  layer 
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7 

ohmic  loss  (electronic)  in  diffusion  layer  =  — 
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(40) 


Eq.  (24)  is  only  valid  for  the  condition  that  vo2at  the 
membrane  I  catalyst-layer  interface  (f  =  0)  is  greater  than 
zero,  which  indicates 

7  <  270{xq2  -  xbN2  [exp(ft7 )  -  1]}  (41) 
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As  the  current  density  exceeds  the  critical  value,  Eq.  (41) 
is  no  longer  satisfied  and  Eq.  (34)  is  not  applicable  either. 
Under  such  a  condition,  another  function  is  proposed  for 
the  dimensionless  oxygen  concentration  profile  within  the 
catalyst  layer,  namely 


f  —  w 

=  XOl  l 


for  w  <  £  <  1  and 
xq2  =0  for  0  <  f  <  w 


(42) 


where  w  denotes  the  dimensionless  position  in  the  catalyst 
layer  where  oxygen  is  completely  consumed  by  the  electro¬ 
chemical  reaction.  Relating  the  transport  rate  of  dissolved 
oxygen  at  the  catalyst-layer  I  gas-diffuser  interface  to  the  elec¬ 
tric  current  density  of  the  cathode  results  in 


2 xs  I 

ZAo2  ' 


Lo 


(43) 


Accordingly, 
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have  the  forms 
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ohmic  loss  (electronic)  of  catalyst  layer  = 


7rs  J 

,xo2l0 

16)01 


/  17 

ohmic  loss  (protonic)  of  catalyst  layer  = - 

Pi  16 
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4.  Approximate  model  II  for  performance  equations 

The  performance  equations  derived  previously  can  be 
further  simplified  if  the  potential  profiles  of  the  solid  cat¬ 
alysts  and  the  ionomer  phase  in  the  catalyst  layer  are  as¬ 
sumed  to  be  piecewise  linear  rather  than  parabolic  as  stated 
in  the  approximate  model  I.  Following  the  postulation 
that  the  electrochemical  reaction  is  lumped  at  the  reaction 
center,  the  ionomer  potential  decreases  linearly  from  the 
membranelcatalyst-layer  interface  to  the  reaction  center  due 
to  an  invariant  protonic  current  density.  In  addition,  the 
ionomer  potential  profile  becomes  flat  after  the  reaction 
center  since  the  protonic  current  density  is  essentially  zero 
in  this  region.  On  the  other  hand,  the  potential  profile  of  the 
solid  catalysts  behaves  oppositely.  It  appears  as  a  horizon¬ 
tal  line  from  the  membranelcatalyst-layer  interface  to  the 
reaction  center  because  of  the  absence  of  electronic  move¬ 
ment  in  this  region  and  decreases  linearly  thereafter  due  to 
a  constant  electronic  current.  Based  on  these  observations, 
the  potential  profiles  for  the  solid  catalysts  and  the  ionomer 
can  be  written  as 


(46) 

(47) 


for  f*  <  £  <  1; 
for  0  <  f 


(52) 


Just  following  the  same  procedures  as  those  without  oxygen 
depletion,  the  cathode  potential  can  be  related  to  the  electric 
current  density  as 
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The  diffusion  overpotential  and  the  electronic  ohmic 
loss  of  the  gas  diffuser  and  the  activation  overpotential  for 
the  electrochemical  reaction  have  the  same  expressions  as 
those  of  Eqs.  (35),  (36)  and  (40).  Other  contributive  terms 


=  --?*  for  f  <  1; 
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=  __L?  for  0  <  £  <  £*  (53) 

Pi 

As  a  result,  the  performance  equation  becomes 
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for  the  condition  without  oxygen  depletion.  Accordingly, 
ohmic  losses  due  to  the  electron  conduction  and  proton 
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migration  are  modified  to 


5.  Results  and  discussion 
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while  other  contributive  terms  in  the  cathode  polarization 
remain  the  same.  As  oxygen  depletion  occurs,  the  perfor¬ 
mance  equation  is  changed  to 
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Thus,  the  ohmic  losses  due  to  electron  conduction  and  proton 
migration  have  to  be  re-evaluated  as 


ohmic  loss  (electronic)  of  catalyst  layer  = 


-V ° 

2/3 1 


(58) 


The  one-dimensional  model  and  the  performance  equa¬ 
tions  formulated  in  the  previous  sections  have  been  used  to 
construct  the  polarization  curves  for  PEFC  cathodes  in  order 
to  investigate  the  accuracy  of  the  performance  equations.  The 
parameter  values  used  for  the  base  case  are  listed  in  Table  1 . 
In  fact,  due  to  the  variations  in  the  techniques  and  in  the 
properties  of  materials  used  for  fabricating  a  PEFC  cathode, 
the  values  of  its  physical  and  chemical  properties  may  differ 
by  some  orders  of  magnitude.  As  mentioned  previously,  the 
major  objective  of  the  present  study  is  to  investigate  the  ac¬ 
curacy  of  the  performance  equations  by  taking  the  numerical 
solutions  of  the  one-dimensional  model  as  exact  ones  for  an 
extended  parameter  range.  In  order  to  do  this,  computations 
were  carried  out  by  adjusting  the  magnitude  of  a  particular 
parameter  in  Table  1,  while  other  parameter  values  remain 
fixed.  The  loci  of  the  performance  equations  as  well  as  the 
solutions  of  the  one-dimensional  model  for  the  base  case 
are  shown  in  Fig.  1.  There  is  rather  good  agreement  among 
these  three  curves  since  they  cannot  be  distinguished  except 
the  region  near  the  limiting  currents.  The  differences  in  the 
limiting  currents  calculated  based  on  the  performance  equa¬ 
tions  and  the  one-dimensional  model  are  caused  partially  by 
the  approximation  (Eq.  (32))  that  is  used  to  evaluate  the  oxy¬ 
gen  concentration  at  the  gas-diffuserlcatalyst-layer  interface, 


ohmic  loss  (protonic)  of  catalyst  layer  =  — - -  — - — 

(59) 

The  performance  equations  derived  based  on  the  approx¬ 
imate  model  II  appear  to  have  simpler  expressions  for  the 
ohmic  losses  pertinent  to  proton  migration  and  electron 
conduction  in  the  catalyst  layer  than  those  based  on  the 
approximate  model  I,  particularly  for  the  condition  that 
oxygen  depletion  does  not  yet  occur.  At  low  current  density, 
however,  such  differences  would  be  of  little  significance 
compared  with  other  contributive  terms  of  polarization.  At 
high  current  densities,  most  dissolved  oxygen  in  the  cat¬ 
alyst  layer  is  consumed  within  a  narrow  region  close  to 
the  catalyst-layerlgas-diffuser  interface.  Consequently,  the 
potential  loss  caused  by  proton  migration  becomes  com¬ 
pletely  dominant  in  the  ohmic  losses  of  the  catalyst  layer 
due  to  much  lower  values  of  protonic  conductivity  and 
much  longer  travel  distance  for  protons.  As  a  result,  the 
differences  between  the  performance  equations  formulated 
according  to  these  two  approximate  models  are  expected 
to  be  negligible.  These  observations  imply  that  the  perfor¬ 
mance  equations  of  the  approximate  model  II  are  derived 
using  simpler  forms  of  potential  profiles  for  the  ionomer 
phase  and  the  solid  catalysts,  but  most  likely  with  no  ap¬ 
preciable  degradation  of  accuracy.  Detailed  comparisons 
based  on  computations  will  be  delineated  in  the  following 
sections. 


Table  1 


Values  of  model  parameters  for  base  case 


Effective  ionic  conductivity  of  ionomer, 

1.06  x 

10“2 

fcjf  (mho  cm-1) 

Effective  electric  conductivity  of  catalyst, 

0.354 

kff  (mho  cm-1) 

Effective  electric  conductivity  of  gas 

0.5 

diffuser,  k^f  (mho  cm-1) 

Effective  diffusivity  of  dissolved  oxygen  in 

4.2  x 

10-7 

the  catalyst  layer,  D 

Effective  gas-pair  diffusivity, 

6.97  x 

10“3 

£>02— n2  (=  £L5£>o2-n2)  (cm2  s-1) 

Effective  gas-pair  diffusivity, 

9.67  x 

10“3 

£>N2f-w  (=  £L5£>n2-w)  (cm2  s-1) 

Effective  gas-pair  diffusivity, 

9.48  x 

10“3 

£>o2-w  (=  £L5£>o2-w)  (cm2  s-1) 

Gas-diffuser  thickness,  d  (pan) 

300 

Catalyst-layer  thickness,  8  (|xm) 

10 

Effective  porosity  of  gas  diffuser,  s 

0.25 

Product  of  platinum  surface  area  and 

4.1  x 

10“4 

reference  exchange  current,  Ai0,ref 

(A  cm-3) 

Pressure,  P  (atm) 

5 

Bulk  oxygen  molar  fraction,  Xq2 

0.190 

Bulk  nitrogen  molar  fraction,  jcjjj2 

0.716 

Molar  fraction  of  water  vapor,  vw 

0.094 

Temperature,  T  (°C) 

80 

Henry’s  constant,  Hq>2  (atm cm3  mol-1) 

2.0  x 

105 

Equilibrium  potential,  V0  (V) 

1.2 

Cathodic  transfer  coefficient,  ac 

0.5 

Anodic  transfer  coefficient,  aa 

0.5 
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I,  A/cm2 

Fig.  1.  Polarization  curves  for  base  case.  (■)  One-dimensional  model; 
(A)  approximate  model  I;  (•)  approximate  model  II. 

and  partially  by  the  numerical  errors  of  the  one-dimensional 
model  due  to  extremely  high  oxygen  concentration  gradi¬ 
ents  developed  while  the  current  density  approaches  its  lim¬ 
iting  value.  The  individual  contributive  terms  of  the  overall 
cathode  polarization  for  the  base  case  shown  in  Fig.  1  are 
displayed  in  Fig.  2,  which  are  calculated  using  the  perfor¬ 
mance  equations  derived  from  the  approximate  model  II.  It 
is  seen  that  the  ohmic  potential  loss  pertinent  to  electron 
motion  in  the  catalyst  layer  is  negligible  for  the  entire  range 
of  current  density  because  of  the  high  electric  conductivity 
of  the  solid  catalysts  and  the  small  thickness  of  the  layer. 
The  activation  overpotential  for  oxygen  reduction  increases 
sharply  at  extremely  low  current  densities.  As  the  current 
density  is  further  increased,  the  activation  overpotential  in¬ 
creases  in  a  more  moderate  fashion  and  other  contributive 
terms  for  the  cathode  potential  losses  gain  more  importance. 
The  potential  losses  for  proton  migration  in  the  ionomer 


Fig.  2.  Polarization  curves  of  individual  contributive  terms  for  the  base 
case  calculated  using  approximate  model  II.  (■)  Ohmic  potential  loss  of 
catalyst;  (O)  activation  overpotential  for  oxygen  reduction  reaction;  (♦) 
diffusion  overpotential  of  catalyst  layer;  (□)  diffusion  overpotential  of 
gas  diffuser;  (•)  ohmic  potential  loss  of  ionomer  in  catalyst  layer;  (A) 
ohmic  potential  loss  of  gas  diffuser. 


Fig.  3.  Polarization  curves  for  cathodes  with  different  catalyst-layer  thick¬ 
nesses.  (A)  5  =  5  [xm;  (B)  8  =  20  [xm;  (■)  one-dimensional  model;  (A) 
approximate  model  I;  (•)  approximate  model  II. 

phase  of  the  catalysts  layer  and  electron  conduction  in  the 
gas  diffuser  are  quite  similar;  both  appear  linear  but  with 
different  slopes.  The  diffusion  overpotential  in  the  catalyst 
layer  exhibits  similar  characteristics  to  those  of  the  activa¬ 
tion  overpotential  at  low  and  intermediate  densities;  never¬ 
theless,  it  increases  sharply  to  infinity  as  the  limiting  current 
is  approached.  Such  a  phenomenon  is  also  observed  for  the 
diffusion  overpotential  for  the  gas  diffuser  due  to  the  vanish¬ 
ing  of  oxygen  concentration  at  the  catalyst-layerlgas-diffuser 
interface. 

The  polarization  curves  for  two  cathodes  with  different 
thickness  of  the  catalyst  layer,  i.e.  5  and  20  (Jim,  are  given 
in  Fig.  3.  It  is  seen  that  the  results  calculated  based  on  the 
performance  equations  of  the  approximate  models  I  and  II 
and  the  one-dimensional  model  are  well  matched.  At  low 
current  densities,  oxygen  is  able  to  reach  every  point  of  the 
catalyst  layer.  Thus,  a  thicker  catalyst  layer  indicates  that 
larger  space  is  available  for  the  oxygen  reduction  reaction. 
As  a  result,  a  lower  activation  overpotential  is  required  to 
generate  the  same  electric  current.  In  other  respects,  higher 
values  of  ohmic  potential  losses  for  proton  migration  is  ex¬ 
pected  for  a  thicker  catalyst  layer  due  to  the  longer  travel  dis¬ 
tance.  Moreover,  a  higher  diffusion  overpotential  is  formed 
by  a  lower  average  value  of  oxygen  concentration  over  the 
catalyst  layer.  Recall  that  potential  loss  caused  by  electron 
conduction  in  the  catalyst  layer  is  negligible.  Computational 
results  reveal  that  the  first  factor  is  more  pronounced  than 
the  other  two  at  low  current  densities,  which  thus  results  in  a 
higher  potential  for  the  cathode  with  a  thicker  catalyst  layer. 
At  intermediate  and  high  current  densities,  however,  oxygen 
depletion  takes  place  within  the  catalyst  layer;  therefore,  the 
influence  of  the  first  factor  diminishes  and  those  of  the  other 
two  become  dominant.  As  a  consequence,  a  higher  poten¬ 
tial  is  found  for  the  cathode  with  a  thinner  catalyst  layer.  It 
should  be  noted  that  the  diffusion  overpotential  of  the  dif¬ 
fuser  does  not  change  with  the  thickness  of  the  catalyst  layer 
as  indicated  by  Eq.  (35). 


34 


H.-K.  Hsuen/  Journal  of  Power  Sources  123  (2003)  26-36 


Fig.  4.  Polarization  curves  for  cathodes  with  different  product  values  of 
Az'o,ref-  (A)  Aiojef  =  4.1  x  10-3  A  cm-3;  (B)  Aioje f  =  4.1  x  10-5  A  cm-3; 
(■)  one-dimensional  model;  (A)  approximate  model  I;  (•)  approximate 
model  II. 

The  polarization  curves  calculated  for  two  cathodes  with 
different  product  values  of  A/0,ref>  viz.  4.1  x  10~3  and  4.1  x 
10~5  A  cm-3,  are  displayed  in  Fig.  4.  Again,  the  polarization 
curves  obtained  using  the  approximate  models  I  and  II  agree 
well  with  those  from  the  one-dimensional  model.  Among 
all  the  parameters  considered  in  the  present  models,  A/0,ref 
is  the  only  one  that  characterizes  the  activity  of  the  cathode 
catalysts.  If  /0,ref  remains  unchanged,  a  larger  product  value 
of  A/0?re f  implies  that  the  cathode  is  composed  of  catalysts 
with  a  higher  platinum  dispersion,  which  provide  a  larger 
active  surface  area  to  accommodate  oxygen  reduction.  For 
a  cathode  with  a  higher  value  of  Aio  re f,  a  lower  activation 
overpotential  is  required  to  generate  an  identical  electric 
current;  thus,  a  higher  cathode  potential  is  expected.  Such 
an  effect  can  be  quantitatively  estimated  by  Eq.  (36).  As 
shown,  a  difference  of  0.07  V  in  the  cathode  potential  results 
from  the  two  different  A/0,ref  values  used  in  calculating  the 
results  presented  in  Fig.  4. 

Investigation  of  the  discrepancies  among  the  polarization 
curves  from  the  approximate  models  I  and  II  and  from  the 
one-dimensional  model  was  carried  out  for  cathodes  with 
two  different  values  of  effective  diffusivity  of  oxygen  in  the 
catalyst  layer  (Z)q^).  The  results  are  presented  in  Fig.  5. 
A  fairly  good  agreement  is  found  between  the  three  mod¬ 
els  as  in  the  previous  cases.  The  parameter provides  a 
quantified  measure  of  the  transport  ability  of  oxygen  un¬ 
der  the  multiphase  structure  of  the  catalyst  layer.  A  higher 
value  of  indicates  that  oxygen  is  able  to  reach  a  deeper 
portion  of  the  catalyst  layer  and,  consequently,  a  lower  dif¬ 
fusion  overpotential  in  the  catalyst  layer  results.  This,  in 
turn,  gives  rise  to  a  higher  cathode  potential.  As  shown,  the 
difference  in  limiting  current  density  calculated  using  the 
one-dimensional  model  and  the  approximate  ones  for  the 
cathode  with  Z)q^  =  2.1  x  10-7  cm2  s_1  is  much  lager  than 

that  with  =  2. 1  x  10-5  cm2  s_1 .  This  is  attributed  to  the 
fact  that  a  steeper  oxygen  concentration  profile  is  formed  for 


Fig.  5.  Polarization  curves  for  cathodes  with  different  values  of  effective 
oxygen  diffusivity  in  catalyst  layer.  (A)  D =  2.1  x  10  5  cm2  s  {;  (B) 

DqJ  =  2.1  x  10~7  cm2  s-1;  (M)  one-dimensional  model;  (A)  approximate 
model  I;  (•)  approximate  model  II. 

the  catalyst  layer  with  a  lower  Z)q^,  which  therefore  leads  to 
larger  numerical  errors  in  the  region  near  the  limiting  cur¬ 
rent.  Nevertheless,  the  present  performance  equations  yield 
the  same  values  of  limiting  current  for  theses  two  cases. 

Polarization  curves  are  constructed  in  Fig.  6  for  two 
cathodes  with  different  values  of  effective  ionomer  con¬ 
ductivity  in  the  catalyst  layer  by  using  the  three  models 
formulated  in  the  present  study.  The  curves  from  the  three 
models  almost  coincide.  In  addition,  the  cathode  with  = 
0.035  mho  cm-1  always  exhibits  a  higher  cathode  potential 
than  that  with  =  0.0035  mho  cm-1  for  the  same  current 
density.  The  computational  results  also  revealed  that  the 
effect  of  on  the  cathode  potential  is  negligible  at  low 
current  densities  and  gains  more  importance  as  the  current 
density  is  increased.  Although  the  expression  of  ohmic 
potential  loss  due  to  proton  migration  involves  exponential 


Fig.  6.  Polarization  curves  for  cathodes  with  different  values  of  effective 
proton  conductivity  of  ionomer  phase  in  catalyst  layer.  (A)  =  0.035 

mho  cm-1;  (B)  =  0.0035  mho  cm-1;  (I)  one-dimensional  model; 

(A)  approximate  model  I;  (•)  approximate  model  II. 
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and  fractional  functions  of  current  density,  it  behaves  in  a 
linear  fashion  as  illustrated  by  the  results  shown  in  Fig.  2. 

The  diffusion  layer  is  a  very  important  component  in 
determining  the  cathode  performance  at  high  current  densi¬ 
ties.  Under  operation,  the  effective  porosity  of  the  cathode 
diffuser  changes  with  the  amount  of  liquid  water  that  ac¬ 
cumulates  in  the  layer.  It  depends  on  the  physical  and 
chemical  properties  of  the  layer,  the  degree  of  external  hu¬ 
midification,  and  the  current  density.  Some  research  efforts 
have  been  made  to  uncover  such  dependencies  [7,10].  Since 
the  emphasis  of  the  present  work  is  placed  on  investigat¬ 
ing  the  accuracy  of  the  developed  performance  equations, 
a  constant  value  of  effective  porosity  is  adopted  for  the 
diffuser  in  calculating  the  polarization  curves.  Once  an  ex¬ 
plicit,  algebraic  expression  that  describes  the  reliance  of  the 
effective  porosity  of  diffuser  on  current  density  and  other 
parameters  is  available  [25],  such  a  relation  can  be  directly 
incorporated  into  the  present  performance  equations.  The 
polarization  curves  for  the  cathodes  with  two  different  s 
values  are  presented  in  Fig.  7,  which  have  also  been  calcu¬ 
lated  using  the  three  models  developed  in  the  present  work. 
As  illustrated,  the  polarization  curves  obtained  from  the 
three  models  are  almost  identical.  At  low  and  intermediate 
current  densities,  where  the  diffusion  limitations  of  oxygen 
within  the  diffuser  are  not  significant,  the  potentials  of  these 
two  cathodes  are  quite  close.  Nevertheless,  at  high  current 
densities,  where  the  diffusion  limitations  of  oxygen  within 
the  diffuser  play  a  more  important  role  in  the  cathode  per¬ 
formance,  the  cathode  with  s  =  0.1  yields  a  lower  potential 
than  that  with  s  =  0.4.  As  revealed  by  Eq.  (32)  and  Table  1, 
variations  of  s  lead  to  different  values  of  oxygen  concentra¬ 
tion  at  the  gas-diffuserlcatalyst-layer  interface.  Thus,  except 
the  activation  overpotential  for  oxygen  reduction,  all  other 
contributive  terms  to  the  overall  cathode  overpotential  will, 
more  or  less,  be  influenced  by  s.  Among  these  terms,  the 
diffusion  overpotential  of  the  diffuser  is  the  one  that  is  most 


Fig.  7.  Polarization  curves  for  cathodes  with  different  effective  porosity 
of  gas  diffuser.  (A)  s  =  0.4;  (B)  s  =  0.1;  (■)  one-dimensional  model; 
(A)  approximate  model  I;  (•)  approximate  model  II. 


effected.  It  is  also  noted  that  the  limiting  current  density 
of  the  cathode  with  s  =  0.1  is  around  0.85  A  cm-2  and  the 
value  for  the  cathode  with  s  =  0.4  is  much  higher,  which 
is  not  covered  by  the  data  in  Fig.  7. 

6.  Summary  and  conclusions 

In  the  present  work,  a  mechanistic  approach  has  been 
presented  to  derive  the  performance  equations  for  PEFC 
cathodes.  In  this  approach,  the  oxygen  concentration  profile 
in  the  catalyst  layer  is  approximated  by  a  parabolic  poly¬ 
nomial  or  a  piecewise  parabolic  relationship,  as  determined 
by  the  occurrence  of  oxygen  depletion.  Two  different  func¬ 
tions,  namely,  a  parabolic  and  a  piecewise  linear  one,  have 
been  suggested  to  approximate  the  profiles  of  the  ionomer 
potential  and  catalyst  potential  in  the  catalyst  layer.  It  has 
been  shown  that  the  piecewise  linear  function  gives  simpler 
expressions  for  the  potential  losses  due  to  electron  con¬ 
duction  and  proton  migration  without  any  penalty  in  the 
accuracy  of  the  cathode  potential.  The  current  density  is 
estimated  though  lumping  the  rate  of  oxygen  reduction  at 
the  reaction  center.  Stefan-Maxwell  multi-component  dif¬ 
fusion  equations  are  employed  for  the  description  of  mass 
transport  of  gaseous  species  within  the  diffuser.  Investiga¬ 
tions  on  the  accuracy  of  the  present  performance  equations 
have  been  carried  out  for  a  wide  range  of  parameter  values. 
Computational  results  reveal  that  the  proposed  equations  are 
capable  of  providing  accurate  predications  for  an  extended 
operation  domain,  namely,  from  a  point  near  the  equilib¬ 
rium  potential  to  the  condition  where  a  limiting  current 
occurs. 
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